Example 3: Curvilinear scalings

<- back to main page

4.3. Example 3: Curvilinear scalings#

Similar to before we consider a two-dimensional example, namely \(\Omega:=\mathbb R^2\setminus B_{R}(\mathbf x_0)\), some center \(\mathbf x_0\in B_1(0)\) and radius \(R\). This time we allow an arbitrary boundary \(\Gamma\) which is convex such that \(B_R(\mathbf x_0)\) is contained in the interior of \(\Gamma\).

We choose curvilinear coordinates given by

\[\mathbf x(\xi,\hat x):=\hat x+\xi\mathbf n(\hat x)\]

(cf. Section 7.1.3).

Again, we create the according mesh and set some parameters.

from ngsolve import *
from netgen.geom2d import *
from nonlin_arnoldis import *
from ngsolve.webgui import Draw
from numpy import array,sqrt,loadtxt
from matplotlib.pyplot import plot,show,xlim,ylim,legend

N = 25         #infinite elements
maxh = 0.1     #mesh-size
sigma = 0.1+0.3j   #complex scaling paramter
order = 5      #fem order
shift = 4-0.5j      #shift for Arnoldi algorithm
center = (0.2,0)    #center of inner circle
R = 0.5            #radius of inner circle

#create geometry
def MakeEggGeometry(Rx,Ryt,Ryb):
    geo = SplineGeometry()
    pts = [(x,y) for y in [-Ryb,0,Ryt] for x in [-Rx,0,Rx]]
    for pt in pts:
        geo.AppendPoint(*pt)

    inds = [(1,2,5),(5,8,7),(7,6,3),(3,0,1)]
    for i in inds:
        geo.Append(['spline3',*i],bc='outer',leftdomain=1,rightdomain=0)
    geo.AddCircle(center,R,leftdomain=0,rightdomain=1,bc='inner')
    return geo

geo = MakeEggGeometry(0.8,0.7,1.2)


#create mesh
mesh = Mesh(geo.GenerateMesh(maxh=maxh))
mesh.Curve(2*order)
Draw(mesh)
BaseWebGuiScene

The weak formulation after a re-scaling of the solution and testfunction in the exterior is given by (7.12).

Contrary to polar coordinates and star-shaped coordinatesthis time we also need the normal vecter \(\mathbf n(\hat x)\) and the curvature \(\kappa(\hat x)\) and the tangential vector \(\mathbf t(\hat x)\).

n = specialcf.normal(2)
kappa = Trace(grad(n))
t = CoefficientFunction((n[1],-n[0]))

For this first simple example we again choose a frequency independent scaling $\(\sigma(\omega):=\sigma_0\in\mathbb C,\)$ (cf. (7.14a)).

We start by creating the large finite element space for implementing our tensor-product method. Contrary to before we need an additional surface space, since (7.12) contains an additional variable to be able to seperate the bilinear forms into a radial and a surface part (cf. Remark 7.1 and Section 7.2).

Gamma = mesh.Boundaries('outer')

fes_int = H1(mesh,order=order,complex=True)
fes_surf = H1(mesh,order=order,complex=True,definedon=Gamma)
fes_surf_l2 = SurfaceL2(mesh,order=order-1,complex=True,definedon=Gamma)

fes_p = ProductSpace(fes_int,*( N*[fes_surf]) )
fes_u = ProductSpace(*( (N+1)*[fes_surf_l2]))

fes = ProductSpace(fes_p,fes_u)

We import the necessary infinte element matrices and prepare the radial matrices:

from infinite_elements import *

ie_mass,ie_laplace,_,ie_mass_x,_,ie_laplace_x,_,_ = ie_matrices(N)

Now we can assemble our bilinear forms.

ds_g = ds(definedon=Gamma)
(p,u),(q,v) = fes.TnT()
p_int,q_int = p[0],q[0]
S = BilinearForm(
    grad(p_int)*grad(q_int)*dx
    
    +sum(1/sigma*ie_laplace[i,j]*p[j]*q[i]*ds_g
       for i in range(N+1) for j in range(N+1) if abs(ie_laplace[i,j])>0)
    
    +sum(ie_laplace_x[i,j]*kappa*p[i]*q[j]*ds_g
       for i in range(N+1) for j in range(N+1) if abs(ie_laplace_x[i,j])>0)
    
    +sum(sigma*ie_mass[i,j]*(u[i].Trace()*(t*q[j].Trace().Deriv()) + (t*p[i].Trace().Deriv())*v[j].Trace())*ds_g
       for i in range(N+1) for j in range(N+1) if abs(ie_mass[i,j])>0)
    
    -sum(sigma*ie_mass[i,j]*u[i].Trace()*v[j].Trace()*ds_g
       for i in range(N+1) for j in range(N+1) if abs(ie_mass[i,j])>0)
    
    -sum(sigma**2*ie_mass_x[i,j]*kappa*u[i].Trace()*v[j].Trace()*ds_g
       for i in range(N+1) for j in range(N+1) if abs(ie_mass_x[i,j])>0)
    
    ,symmetric=True).Assemble()

M = BilinearForm(
    -p_int*q_int*dx
    -sum(sigma*ie_mass[i,j]*p[j]*q[i]*ds_g
       for i in range(N+1) for j in range(N+1) if abs(ie_mass[i,j])>0)
    
    -sum(sigma**2*ie_mass_x[i,j]*kappa*p[j]*q[i]*ds_g
       for i in range(N+1) for j in range(N+1) if abs(ie_mass_x[i,j])>0)
    
    ,symmetric=True).Assemble()

Finally, we solve the resulting eigenvalue problem.

gf = GridFunction(fes,multidim=100)

#lam = sqrt(array(ArnoldiSolver(S.mat,M.mat,freedofs=fes.FreeDofs(),vecs=gf.vecs,shift=shift**2)))
lam = sqrt(array(PolyArnoldiSolver([S.mat,M.mat],shift**2,200,nevals=80,vecs=gf.vecs,inversetype='sparsecholesky',freedofs=fes.FreeDofs())))

plot(lam.real,lam.imag,'x',label='resonances')
plot([0,5*(1/sigma).real],[0,5*(1/sigma).imag],label='essential spectrum')

#load reference resonances from file
loaded=loadtxt('dhankel_1_zeros.out')
ref=(loaded[:,0]+1j*loaded[:,1])/R

plot(ref.real,ref.imag,'.k',label='reference')

xlim((0,10))
ylim((-5,0))
legend()
Draw(gf.components[0].components[0])
1/200
2/200
3/200
4/200
5/200
6/200
7/200
8/200
9/200
10/200
11/200
12/200
13/200
14/200
15/200
16/200
17/200
18/200
19/200
20/200
21/200
22/200
23/200
24/200
25/200
26/200
27/200
28/200
29/200
30/200
31/200
32/200
33/200
34/200
35/200
36/200
37/200
38/200
39/200
40/200
41/200
42/200
43/200
44/200
45/200
46/200
47/200
48/200
49/200
50/200
51/200
52/200
53/200
54/200
55/200
56/200
57/200
58/200
59/200
60/200
61/200
62/200
63/200
64/200
65/200
66/200
67/200
68/200
69/200
70/200
71/200
72/200
73/200
74/200
75/200
76/200
77/200
78/200
79/200
80/200
81/200
82/200
83/200
84/200
85/200
86/200
87/200
88/200
89/200
90/200
91/200
92/200
93/200
94/200
95/200
96/200
97/200
98/200
99/200
100/200
101/200
102/200
103/200
104/200
105/200
106/200
107/200
108/200
109/200
110/200
111/200
112/200
113/200
114/200
115/200
116/200
117/200
118/200
119/200
120/200
121/200
122/200
123/200
124/200
125/200
126/200
127/200
128/200
129/200
130/200
131/200
132/200
133/200
134/200
135/200
136/200
137/200
138/200
139/200
140/200
141/200
142/200
143/200
144/200
145/200
146/200
147/200
148/200
149/200
150/200
151/200
152/200
153/200
154/200
155/200
156/200
157/200
158/200
159/200
160/200
161/200
162/200
163/200
164/200
165/200
166/200
167/200
168/200
169/200
170/200
171/200
172/200
173/200
174/200
175/200
176/200
177/200
178/200
179/200
180/200
181/200
182/200
183/200
184/200
185/200
186/200
187/200
188/200
189/200
190/200
191/200
192/200
193/200
194/200
195/200
196/200
197/200
198/200
199/200
200/200
BaseWebGuiScene
_images/32fa6752b9cefa1351b0287079592375bdad5e9ae732670d08d130df04bd541c.png

<- back to main page